Two-dimensional dynamics of expansion of a degenerate Bose gas 
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Expansion of a degenerate Bose gas released from a pancake-like trap is numerically simulated 
under assumption of separation of motion in the plane of the loose initial trapping and the motion 
in the direction of the initial tight trapping. The initial conditions for the phase fluctuations are 
generated using the extension to the two-dimensional case of the description of the phase noise by 
the Ornstein-Uhlenbeck stochastic process. The numerical simulations, taking into account both 
the finite size of the two-dimensional system and the atomic interactions, which cannot be neglected 
on the early stage of expansion, did not reproduce the scaling law for the peaks in the density 
fluctuation spectra experimentally observed by Choi, Seo, Kwon, and Shin [Phys. Rev. Lett. 109, 
125301 (2012)]. The latter experimental results results may thus require an explanation beyond our 
current assumptions. 

PACS numbers: 67.85.-d,03.75.Lm,03.75.Hh 
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Correlations in degenerate Bose gases with repul- 
sive interactions manifest themselves most apparently 
via phase coherence. Even low-dimensional degenerate 
bosonic systems (quasicondensates) demonstrate phase 
coherence over finite distances. After release from the 
trap, phase fluctuations are converted, in the course of 
free expansion, into density-density correlations (density 
ripples). Theory for expansion of both one-dimensional 
(ID) and two-dimensional (2D) trapping geometries has 
been developed [l| and experimentally proven for ID 
trapped gases [2| (similar behavior has been previously 
found in expanding clouds, which were in the three- 
dimensional Bose-Einstein condensate rather in the ID 
regime inside a very elongated trap Q). 

In contrast to ID Bose gases, where correlations decay 
exponentially and their characteristic length At deter- 
mines the typical wavelength scale of the density ripples 
and the timescale mX^/h of their emergence (m being 
the mass of the atom) , correlations in degenerate 2D Bose 
gases decay with the distance according to a power law, 
if the temperature is below the point of the Beresinskii- 
Kosterlitz-Thouless (BKT) transition Such a power- 
law decay renders no specific correlation length, the spec- 
trum of the density fluctuations of a 2D gas released from 
a trap evolves in a self-similar way at asymptotically large 
expansion times t c . In particular, after averaging the 
power spectrum of density fluctuations 
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where 5n(p) is the local density fluctuation, over the di- 
rections of the 2D wave vector q, the nth peak position 
q n , n = 1, 2, 3, . . . , is expected to satisfy the relation [l[ 



% 2 ! t e /(27rm) = n - 1/2. 



(2) 



Recently, Choi et at 5] reported their experimental 
results on probing fluctuations in a 2D gas of bosonic 
( 23 Na) atoms in free expansion after a sudden release 
from a pancake-like optical trap, the trapping frequencies 



being equal to (u> x , uj v , uj z ) = 2n x (3.0, 3.9, 370) Hz. 
The chemical potential was about (iq — 2irhx 260 Hz, i.e., 
less than the spacing between the discrete levels of the 
trapping Hamiltonian in the tighgtly confined direction. 
The coherent fraction of atomic ensemble decreased from 
0.78 to 0.12 as the temperaure T grew from 20 nK up to 
the BKT transition temperature Tbkt ~ 70 nK. 

The experimental results @ turned out to be quite 
surprising. Instead of recovering the scaling ([5]), they 
brought about another dependence, 



hq^ l t c /(2TTm) = am? 
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with 0.2 < a < 0.45 and 0.7 < 7 < 1 for expansion times 
in the range 10 ms < t c < 25 ms. In other words, the 
spectral peaks of the density fluctuations were not only 
shifted with respect to Eq. but the spacing between 
the ajacent peaks was significantly smaller than Eq. @ 
predicts. Choi et al. ruled out the effects of the cloud fi- 
nite size and suggested that such an unexpected behavior 
might be intrinsic to the expansion dynamics @. 

To check this assumption, we performed numerical 
simulations of 2D expansion of a degenerate gas. Since 
the temperature and, hence, the coherent fraction do not 
affect significantly the peak positions q n [f|, we find it 
safe to neglect thermal population of the excited lev- 
els of motion in the tightly trapped (z) direction and 
model the system's dynamics by the Gross-Pitaevskii 
equation (CPE), taking into account the interaction ef- 
fects, which are important at the initial stage of expan- 
sion. Furthemore, we separate the motion in z-direction 
and in the [x, y)-plane. This would be impossible if the 
Beliaev and Landau damping of elementary excitations 
in expanding but still dense enough atomic cloud were 
giving rise to significant population of modes with non- 
zero z-components of kinetic momentum. However, the 
estimation of the ratio of the damping rate of an ele- 
mentary excitation to its frequency yields @ the value 

~ aV 2 1 \n\^ A 2 h ) ~ 10~ 2 , where a s is the atomic s-wave 
scattering length, 7130 is the three-dimensional atomic 
number density, and Ath is the de Broglie wave length of 
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atoms at temperature T. Even for elementary excitations 
with the energy close to the chemical potential the typi- 
cal relaxation time (and, hence, the characteristic time of 
scattering into z-direction) in a three-dimensional degen- 
erate gas with the bulk density equal to the peak density 
in our 2D system is of the order of 100 ms, which is much 
longer than the typical expansion time. Therefore we can 
separate the motion in different directions. In z direction, 
we have a free expansion of a wave packet, which is ini- 
tially, at t — 0, the wave function of the ground state 
of the harmonic trapping potential with the fundamen- 
tal frequency u z . Then the motion in the {x, y) plane is 
described by the GPE 
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The fast (but still not infinitely fast) evolution of 
the quasicondensate profile in z-direction affects the 
slow motion in the perpendicular plane) only through 
the time-dependent nonlinearity in Eq. The 
time-dependent effective 2D coupling constant fi(t) 
takes into account the growth of the transverse pro- 
file width oc yl + (u> z t) 2 . The normalization condition 
max|*(x,y,t = 0)| = \V(x = 0, y = 0, t = 0)| = 1 fixes 
the prefactor, and we obtain 



/i(f) 



Mo 



(5) 



Initial shape of \^\ 2 is the inverted-parabolic Thomas- 
Fermi profile corresponding to the chemical potential /io 
of the trapped gas of sodium atoms. The phase <p = 
arg 9 at t = can be generated using the generalization 
of the stochastic Ornstein-Uhlenbeck process, previously 
used to model random phase distributions of phases in 
ID quasicondensates pi Iff. 

Assume that the temperature T of 2D Bose gas is below 
the BKT transition, phase and density fluctuations are 
small and, therefore, the Hamiltonian can be linearized: 



H = — 

2m 



(Vi Sn) 
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Here n2r> — J dz is the 2D number density, g is the 
dimensionless coupling strength (g « 0.013 in the experi- 
ment (H), V_l is the 2D gradient in the (x, y)-plane. Then 
the partition function of a trapped gas is 



Z= JvfJ V5n exp[-H/(k B T)]. 



(7) 



The density fluctuation variable can be integrated out, 
and we obtain 



Z = const / T><p exp 



2mk^T 



J d 2 pn 2D {p){\7 ± i P y 



(8) 

We can use the thermodynamic expression © to gen- 
erate initial phase distribution for individual realizations 



of our numerical simulations. On a (2M + 1) x (2M + 1) 
square grid Eq. © takes the form 



M M „ r 

Z = const ]J II / d(fil * l v exp { 

l lc =-Mly = -M'' { 
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[ki s - Wx/y-i) 2 + - ¥>/*-nJ 2 ] j- ( 9 ) 

Here the pair of subscripts l x and l y denotes a vari- 
able at the point with coordinates x — l x As and 
y = l y As [note that the grid step As does not ap- 
pear in Eq. ([9]) in our 2D problem]. The density 
profile and its peak value are characterized by dimen- 
sionless quantities wi x i y = n2v>{l x As, l y As)/ri2T)(Q, 0) 
and e = ?i 2 n2D(0, Q) / (mkB,T) , respectively. Periodic 
boundary conditions are set, i.e., tpi x m+i = <pi a —m and 

Assume we know the 2M + 1 phases ipi x —i i ■ Then we 
can calculate the phases ipi x i in the next row of the grid 
using the transfer matrix approach [9[. 

We assume that the grid step is small enough and thus 
we can write wi x i y /wi x0 s» wi x -n y /wi x -i . The trans- 
formation [hi 
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K i:i= E V ^r ' 



(10) 



l v =-M 

diagonalizes the non-negative quadratic form 
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l y = -M '*° j=Q 

(11) 

Both Vi x; ji y and Ji x j can found numerically by standard 
methods of linear algebra. Then Eq. ^ factorizes: Z — 
const n 2 =o Zj, where each factor 
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corresponds to a ID evolution (along the x-coordinate) 
of a Aj describable by the Ornstein-Uhlenbeck stochas- 
tic process (which can be easily seen from the anal- 
ogy with the relative phase between two tunnel-coupled 
ID quasicondensates at thermal equilibrium |7[). Here 
ei = ewi q. The formula [111 
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where <; is a (pseudo)random number obeying the Gaus- 
sian distribution with zero mean and unity variance, up- 
dates the variable Aj on the next step. Then we obtain 
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FIG. 1: Solid line: the power spectrum of density fluctuations 
for t c — 11 ms obtained from Eq. Q, i.e., taking into account 
interaction of atoms in an expanding atomic cloud. Parame- 
ters of the trapped degenerate cloud are taken from Ref. [H 
(see details in the text). Dashed line: spectrum of density 
fluctuations obtained for a purely ballistic expansion during 
the same t e . The spectra are averaged over 5 realisations each. 

the initial phase values in the next grid row: 

M 

<Pt.l,= E V ^ly\7^- A ^r ( 14 ) 

ly=-M V l " l V 

Subsequent steps are done using the new local values for 
the density profile. 

Eq. (|13jl. together with chosing the (pseudo)random 
initial values of the Aj variables according to a Gaus- 
sian distribution with zero mean and variance (Af ■ ) = 

l/(2e/ a . yj J; x j), enables us to prepare the initial condi- 
tions for each realization of our numerical modeling [l2j . 

We solved Eq. (J4j) using the split-step pseudospectral 
method for the typical parameters of Ref. [j| and ob- 
tained the power spectra of density fluctuations. As the 
density fluctuation Sn in Eq. (JT|) we took the difference 
between the local densities given at the same expansion 
time t c by two solutions, one containing initial phase fluc- 
tuations and other having constant phase everywhere at 
t = 0. An example of such a spectrum P q , that is P q aver- 
aged over the direction of the wave vector q, is displayed 
in Fig. [TJ The variation of the temperature from 60 nK 
to 20 nK (corresponding to the increase of e from 15 to 
45) changes the height of the peaks of P q only, leaving 
the peak positions q n unmodified, as we can expect from 
the experimental results @. The solution of the GPE g]) 
is also juxtaposed with the density fluctuations spectra 
for purely ballistic expansion [i.e., by setting fj,(t) = 0]. 
In the case of purely ballistic expansion, we found that 
Eq. @ is satisfied with 1% accuracy for n > 2; for the 
first peak position we obtained hqft e / (27rm) ps 0.4 for 
10 ms < t c < 23 ms. 

As we can see from Fig. [TJ interaction of atoms during 
expansion is not negligible and shifts the peak positions 
with respect to the ballistic expansion case. However, 



instead of recovering Eq. (j3]), we found that the peak 
positions are well described by the formula 

ft«*fo/(27rm) = P{n - 1), (15) 

where 

/3«0.9, J?«0.7 (16) 

for the whole range of experimentally relevant expansion 
times (see Fig. [5]). 

Note, that the transversal trapping frequency uj z — 
2ir x 370 Hz of Ref. [j| is approximately equal to the 
radial trapping frequency in the experiment of Ref. Q , 
but the cigar-shaped configuration of the atomic cloud 
resulted in the latter case in a faster (oc t~ 2 instead of 
cx i -1 at expansion times > 1 ms) decrease of the three- 
dimensional density of atoms than in the 2D case. Ra- 
dial trapping frequencies in recent experiments with ID 
ultracold gases on atom chips 0, [1] were by an order of 
magnitude higher, thus further decreasing the time scale, 
on which atomic interactions still play a role. Therefore 
expansion of atomic clouds in the experiments @, H ||| 
was very close to ballistic, in contrast to the case modeled 
in our present work. 

Our numerical results do not reproduce the observa- 
tions by Choi et al. Q. For example, our calculations 
predict the values of the 5th peak position 95, which are 
at relevant expansion times larger by a factor 2-3 than 
the observed ones. In general, our calculations predict 
the linear depenfdence of g„ on the peak number n, in 
contrast to the nonlinear fitting formula ((3|) found in Ref. 

Such a conspicuous disagreement between the theory 
that assumes separation of motion in the former trap 
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FIG. 2: The peak positions for the power spectrum of density 
fluctuations (the peak number n is dimensionless) of a 23 Na 
ultracold gas interacting during the expansion for t e — 11 ms 
(filled circles), 14 ms (filled squares), 17 ms (crosses), 20 ms 
(open squares), and 23 ms (triangles). Straight lines (solid, 
long-dashed, dotted, short-dashed, and dot-dashed, respec- 
tively) display the fitting by Eqs. (fTSlfTrl)) . 
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plane and in the direction of the tight confinement be- 
fore the release of the atom from the trap and the ex- 
perimental results Q calls for explanantion. We do not 
expect the scattering of phonons into z-direction in the 
expanding quasicondensate to be a reason for such a dis- 
crepancy, since the estimations based on the theory of 
Ref. @ show that this process is too slow compared to 
the typical expansion times. Moreover, the influence of 
such a scattering would grow with t e , but our results de- 
viate strongly from the experimental findins at shorter 



times as strongly as at longer ones. 

The fact that the experimental spectra of density fluc- 
tuations are describable by Eq. ^ when no vortices 
are present in the 2D system [j| allow us to rule out 
any vortex-related explanation of this discrepancy, which 
therefore still remains puzzling. 
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